Vortex solitons in an off-resonant Raman medium 
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We investigate existence and linear stability of coupled vortex solitons supported by cascaded four- 
wave mixing in a Raman active medium excited away from the resonance. We present a detailed 
analysis for the two- and three-component vortex solitons and demonstrate the formation of stable 
and unstable vortex solitons, and associated spatio-temporal helical beams, under the conditions of 
the simultaneous frequency and vortex comb generation. 



PACS numbers: 

I. INTRODUCTION 

Optical vortices are point phase singularities of the 
electromagnetic field, with the beam intensity vanishing 
at the singularity and the field phase changing by 27r/ 
along any closed loop around it. / = 0,±1,±2,... is 
known as the orbital angular momentum quantum num- 
ber or vortex charge. In a nonlinear medium vortices can 
propagate undistorted due to a balance between diffrac- 
tion and nonlinearity, and form so-called vortex solitons 
[T]. Nonlinearity can also trigger frequency conversion 
accompanied by the conversion of the charge In par- 
ticular, in the second harmonic generation process, the 
fundamental field carrying a vortex with the charge / is 
converted into the second harmonic field with the charge 
2/ dJ 121 |3l m [5l [6] . Analogous conversion rules have been 
reported for the degenerate four-wave mixing in Kerr- 
like materials [7] and for the three- wave Raman resonant 
process [8^. Multi-component vortex solitons sustained 
by the interaction of the beams with different frequen- 
cies in both quadratic and cubic materials are also well 
known, though under the most typical conditions the fi- 
nite radius vortex solitons break into filaments due to 
azimuthal instabilities [H [21 [4j |5] . 

While the above mentioned experimental and theoret- 
ical research of nonlinear vortex charge conversion has 
focused on cases involving a small number of frequency 
components, typically two or three, the efforts directed 
towards short pulse generation have resulted in the devel- 
opment of techniques leading to the generation of dozens 
of coherent frequency side-bands, by means of cascaded 
four- wave mixing in Raman active gases [9 , 10 . The lat- 
ter technique does not rely on the waveguide or cavity 
geometries to boost nonlinear interaction and is there- 
fore suitable for the simultaneous frequency and vortex 
charge conversion. This idea has been explored by our 
group and we have recently demonstrated simultaneous 
generation of frequency and vortex combs [11] in a Ra- 
man medium excited off-resonance with the two pump 
beams, when one of the two carries a unit vortex and 
the other is vortex free. We have derived the vortex 
conversion rules and demonstrated that the simultane- 



ous frequency and vortex combs are shaped in the form 
of the spatio-temporal helical beams [11]. On the focus- 
ing side of the Raman resonance, the mult i- component 
vortex solitons have been found. 

The aim of this work is to report regular tracing of the 
multi-component vortex solitons in the parameter space 
and to study their linear stability with respect to per- 
turbations. Our analysis shows that the spectrally sym- 
metric soliton solutions centered around the vortex-free 
frequency component are typically unstable, although the 
instability fully develops only after long propagation dis- 
tances. At the same time, the asymmetric solitons, for 
example those where all the generated components are 
the Stokes ones, have a broad stability range. Based 
on the results of the linear stability analysis for 2 and 
3 component solitons, we demonstrate the same general 
tendencies of the soliton dynamics for the case of many 
coupled side-bands. 



II. MODEL 

The dimensionless model describing the evolution 
of the side-bands in an off-resonant ly excited Raman 
medium is [3 |TT] 

id.Er, - ^AEn = PnEn + Q*^n-1 + (1) 

where n = -M + 1, . . . , 0, . . . , A/" (M, A" > 0), and A = 
^ ^y' En are the dimensionless amplitudes of the 
sidebands, such that the total field is given by 

Etot = Y.E^{x,y,z)^''-'-'''-\ (2) 

n 

where Vtn (^o + nUrnod) /^mod, ^mod = cji - cjq IS the 
modulation frequency (i.e. the frequency difference be- 
tween the two driving fields). A" is the number of the 
anti- Stokes components and M — 1 is the number of the 
Stokes components. Taking into account the Eq field, 
we have M -\- N interacting Raman side-bands. The 
physical frequencies and wavenumbers are represented 
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FIG. 1: (color online) . Off-resonant excitations of the Raman 
transition creating either focusing or defocusing nonlineari- 
ties. cOr is the Raman frequency and comod — ^\ — where 
uj\ > ooo are the pump frequencies. 



by the lower case letters ujn and /c^, whilst their di- 
mensionless counterparts by the upper case: and 
Kn. The dimensionless time t is measured in units 
of 1/uJmod^ the propagation coordinate z is in units of 
L, and the transverse coordinates (x^y) are in units of 
yZc/oJo- = {(JQ + nUmod)L/ c are the scaled free 

space wavenumbers. Here, L = {r]huJoJV\b\)~^ character- 
izes the coupling length over which power is transferred 
between neighboring side-bands in the absence of disper- 
sion. T] ^ 376 is the free space impedance, JV is the 
density of molecules and 6 is a coefficient characterizing 
the material dependent coupling between the sidebands 
[9]. The weak frequency dependence of b is neglected for 
simplicity. 

L varies from 1 to a few mm for D2 and H2 gases ^ , so 
that one unit of x corresponds to a few tens of microns. 
Q is the Raman coherence responsible for the coupling 
between the side-bands. Neglecting dissipation due to 
finite linewidth of atomic transition and finite dephasing 
time, in the adiabatic approximation ^ [Hi [T2] 



Q{En) 



sgn{jii)S 



where ja = IcOmod — ^r|/(|^|^o) is the scaled modulus of 
the detuning of the modulation frequency from the Ra- 
man frequency ujr- We also note, that the above result 
is obtained under the assumption of equal Stark shifts of 
molecular levels, which is the case for large detunings [9^. 
While can always be fixed to unity by proper rescaling 
of the field amplitudes, its sign controls the effective type 
of nonlinearity in Eqs. ([T]): positive (negative) /i corre- 
sponds to the focusing (defocusing) nonlinearity [HI [13], 
see Fig. [T] In what follows we consider the case of the 
focusing nonlinearity [13, J^, (/i = 1), which is known 
to support bright soliton solutions [T5j [16] . 

IQI varies from to 1/2 for varying from to 00. 

Therefore nonlinear interaction between harmonics is sat- 
urated at high powers or, equivalently, at small detunings 
l^mod — ^r\' ^ua/Tq are the dimensional amplitudes of 
the harmonics. For D2 and H2 gases /i = 1 corresponds 
to /q ^ O.lGW/cm^, provided {oOmod — ^ IGHz. 
/3n = I3{(j0n) is the propagation constant of the nth har- 



III. 



SOLITON SOLUTIONS: 
FRAMEWORK 



GENERAL 



In this and the next chapter we describe the gen- 
eral framework for finding the stationary soliton solu- 
tions and studying their linear stability. Application of 
these techniques to the cases of two and three compo- 
nents are described in detail in Chapters V and VI. The 
fact that equations ([T]), ([3| are invariant with respect to 
En En exp(z^) and En En exp(mV^), where (j) and ijj 
are arbitrary constants [TTl [12] , implies the conservation 
of the two integrals P = In and R = nin-, where 
In = J J dxdy\En\^^ and suggests the following ansatz 
for the soliton solutions: 

En{x, y, z) = fn{r) exp [ilnO + + K2n)z] . (4) 

Here r and 6 are the polar radius and angle, /n = ^0 + 
^{h—lo) is the vortex charge of the nth harmonic, hii^2 are 
free parameters associated with the above symmetries. 
The choice of Iq and li defines the step, Al = li — Iq, in 
which the vortex charge is changing between the adjacent 
side-bands, /n(^) are real functions obeying 



r dr 



JL f 



[kI + H2n + Pn] fn 
<7(/n-l +/n+l) . (5) 



where q = Q{fn)' The boundary conditions are [4j: 

fn{r) ^ ci°)rl'"l, r^O, (6) 
fn{r) 



(7) 



where cl^'^^ are real constants. Eq. (6) naturally implies 
that the amplitude of a vortex carrying component. In 7^ 
0, is zero at the phase singularity and that the vortex free 
components, = 0, reach some constant value at r = 0. 
For the fields to decay to zero at r ^ 00, one needs to 
select K\^2 to satisfy 



(8) 



simultaneously for all n. Without any loss of general- 
ity /3o can always be set to zero by the rotation of the 
common phase [TT. Thus, fixing /3o = we find that the 
above inequality for n = gives k\ < 0. At the boundary 
points of the above conditions jn tends to zero. Detun- 
ing the K\^2 values away from these boundaries into the 
range where In is increasing eventually leads to the co- 
herence tending to its maximal value q = 1/2. Examples 
of the radial profiles of the vortex solitons are shown in 
Fig. [2] for the asymmetric configuration with only Stokes 
components being excited {N = 0). When the propa- 
gation constants /3n are symmetric around the central 
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into Eqs. ([T]). After linearization we find: 
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(10) 



FIG. 2: (color online). Soliton radial profiles for the config- 
uration with only Stokes components being excited (A^ = 0), 
(3n — 0.005n^ (which corresponds to normal dispersion), 
charge conversion step is unity: In = n. (a) 2 coupled fields 
(M = 2), Ati = -0.25, = -0.4 (black curves), ^2 = -0.7 
(red/gray curves). Solid (dashed) curves correspond to /o 
(/-i); (b) 5 coupled fields (M = 5), ki = -0.25, K2 = -0.7. 



component: (3^ = P-ni Eqs. ([5| are invariant under the 
transformation n — n, K2 —1^2^ In —In- In this 
case solitons in the opposite configuration, with only anti- 
Stokes components being excited (M = 1), have exactly 
the same structure as those shown in Fig. [2] 

In the vortex soliton case the q = 1/2 limit is 
achieved not only through the growing amplitudes, but 
also through the expansion of the rings and flattening 
of their profiles. The soliton existence boundary corre- 
sponding to q = 1/2 can be worked out neglecting the 
X and y dependence of E^. Subsequently one can dis- 
regard the left-hand sides in Eq. ([5|, assume q = 1/2 
and work out constraints on the hzi and values from 
the solvability conditions of the resulting homogeneous 
equations: [a^i + /^2^ + /5n]/n + (/n-i + /n+i)/2 = 0. Ex- 
amples of the existence domains in the {hii^ /T:2)-plane can 
be seen in Figs. [3]; a) and[7l 



IV. LINEAR STABILITY ANALYSIS: GENERAL 
FRAMEWORK 

Stability of the vortex solutions is of course an impor- 
tant problem, since similar solutions in other models are 
known to exhibit strong modulational instability along 
the rings [H |5l [17] . This instability can be suppressed by 
nonlocal nonlinearities [T8[ [19] , and in some cases when 
the higher order nonlinearities (e.g. quintic) are assumed 
to dominate over the lower order ones (e.g. cubic), see, 
e.g., [71 [20]. Our model is particularly interesting be- 
cause, as we will demonstrate below, it allows the ex- 
istence of a sufficiently broad parameter range, where 
stable vortex solitons exist with the local type of non- 
linearity derived from the first principles. The latter is 
true since the nonlinearity in Eq. ([T]) is calculated from 
the Schrodinger equation for a Raman medium driven far 
from the resonance [O] [12 . 

In order to analyze the linear stability we add small 
perturbations to the vortex solitons and substitute 
the following ansatz 

En = [fn{r)^en{r,6,z)]exp[i{Ki^K2n)z^ilnO] . (9) 



where 



n, m : 





= fn-lMm + fn+lPk , 




= fn-lPm + /n+l^m 7 


dq"" 


_ (2/i2 + 5^) fm-1 - S^fm+1 


dfm 


4(/i2+52)^/' 


dq 


_ (2/i2 + S^) fm+l - S^fm-1 


dfm 


4(/i2+52)^/' 


f 1,.. 


.,0,...iVand8 = 5(/n). 



(11) 

(12) 
(13) 

(14) 



Expanding perturbations into azimutal harmonics [4] 
en{r,0,z) = ^ |/i+j(r,z)exp(iJ(9) 

J>0 

+ {h-j{r,z)rexp{-iJ0)} , (15) 

we assume h^j{r,z) = j{r) ex.-p{Xjz) and derive the 
eigenvalue problem: 
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(16) 



where = {gf_M^j , dt-M^j , • • • , , • • • , Qn^jV - 
and B are the (A/' + M) x (A/' + M) matrix operators. 
Elements of B are Bnm and they are defined in Eq. (12), 
and the elements of are 



d I d 

r 



2r2 



2r dr \ dr ^ 

(^n+l,m ~^ ^n— l,m) ^nm -> 



(17) 



where 5n^m is the Kronecker symbol. For a solution fn to 
be linearly unstable there must exist Aj with Re{Xj) > 0. 
Boundary conditions for eigenstates gj are defined in 
a similar way to the boundary conditions for fn (see 
Eqs. (|6|-([7|), but with In being replaced by In ± J. We 
solve the eigenvalue problem in Eq. ( 16 ) numerically, re- 



placing differential operators by the second-order finite 
differences. Note that accurate stability analysis of the 
multi-component solutions is rather complicated. There- 
fore we will reveal basic mechanisms of instabilities of 
coupled vortex solitons by focusing on two- and three- 
component configurations. Then we will demonstrate by 
numerical modeling of Eqs. ([T]), that the instability and 
stabilization mechanisms found in the simplest cases can 
be seen in the multi-component dynamics. 
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V. TWO-COMPONENT VORTEX SOLITONS 

We start with the simplest configuration of two side- 
bands, that is n = 0, 1 (TV = M = 1) in Eqs. This 
apphes e.g. to the opposite circularly polarized driving 
fields £^0 and Ei , when the cascaded generation of Stokes 
and anti-Stokes harmonics is forbidden due to angular 
momentum selection rules [15j. The propagation equa- 
tions in this case are 



id^.lA-Po)Eo = , (18) 



2 ''V 2^/i^^T\WW 



■ (19) 



Eqs. (18, 19) explicitly express a known fact that the 
fields interacting via the Raman nonlinearity do not have 
nonlinear self-action. This property does not depend on 
the number of interacting components. 

Bright (vortex free) spatial solitons in the two- 
component Raman model have been studied in [151 Ull 
and the associated self-focusing effects have been ob- 
served in [13 , 14 . Also, there are closely related recent re- 
sults on spatial solitons in Raman active liquids [lU [22] . 
The papers [23 have analyzed the two-component tem- 
poral Raman solitons existing in the presence of group 
velocity dispersion, i.e. when the transverse Laplacian is 
replaced with the 2nd-order time derivative. The above 
two-component model is also similar to that for the so- 
called holographic solitons [24] . 

The existence conditions for soliton solutions in 
Eqs. ([8| are reduced to the joint inequalities hii < and 
K:2 < —1^1 — pi ^ which define a semi-infinite region in 
(/^i, K:2) bounded by the two rays, see Fig.jsja). Another 
boundary is derived from the q = 1/2 condition and is 
given by: 



1^2 > 1^2^'^ = 



4a^i 



^1 - Pi 



(20) 



Our linear stability analysis demonstrates that the soli- 
ton with lo = and /i = 1 is unstable only inside the suf- 
ficiently narrow range of K:i^2 corresponding to the rela- 
tively small values of see Fig. [3] As soon as q increases 
and the saturation effects become important the solution 
becomes stable. Note that the saturation of the self- 
focusing nonlinearity does not stabilize the vortices in 
the models with the nonlinear self-action effects [4 . It 
suggests that the absence of the self-action plays an im- 
portant role in stabilization of the vortex solitons. In its 
instability range, the vortex soliton is unstable with two 
eigenvalues Xj=i and Xj=2 having positive real parts, see 
Fig. [4] Fixing a^i we numerically find the critical values 
of f<i2^ at which the two instabilities disappear, see circles 
and squares in Fig.jsj^a). Selective numerical runs for the 
cases lo = 0, |/i| > 2 and Iq = 1, |/i| > 1 suggest that 
they are unstable with respect to azimuthal instabilities 
through large parts of their existence domains. 




(b) 




FIG. 3: (color online). (a) Region of existence of cou- 
pled vortex solitons for the case of two fields, Pi = 0.005. 
Straight blue lines correspond to the boundaries of existence 
in Eq. ([s]), red/gray curve - to the dispersion ^2*'^*'* of high- 
intensity constant amplitude waves, Eq. (20). Shaded area 



indicates region of unstable solutions for the configuration 
{lo = 0,/i = 1), open squares and filled circles correspond 
to numerically found instability thresholds for the J = 1 and 
J = 2 unstable perturbations, respectively (see main body 
text for details and Figs. [4] [s] [6|. (b) Soliton power P and 
maximum value of the coherence q versus K2 at fixed value of 
Ki — —0.25. Approaching the boundary K2^"'^\ki) q is satu- 
rated at its maximum value |^| = 0.5, while the norm tends 
to infinity. 




FIG. 4: (color online). Real part of eigenvalues corresponding 
to unstable perturbations versus k,2 for ki = —0.25. 



To reveal the impact of instabilities on the soliton dy- 
namics, we initialize Eqs. ([T]) with numerically found soli- 
ton solutions slightly perturbed along unstable eigenvec- 
tors and perform dynamical simulations. Results are pre- 
sented in Figs. [5] and [6] for the J = 1 and J = 2 unstable 
eigenvectors, respectively. Both perturbations break the 
soliton symmetry and eventually lead to the formation of 
a single or a pair of bright spatial solitons [151 US] • 



VI. THREE-COMPONENT VORTEX SOLITONS 

The addition of the third component makes the in- 
teraction between the Raman side-bands phase-sensitive, 
and the choice of the vortex charges In in any two fields 
defines the charge of the remaining field via the phase- 
matching conditions [11]. Eqs. ([T]) for the three compo- 
nent case with n = 1 — M, 2 — M, 3 — M are: 
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FIG. 5: (color online). Dynamics of the unstable soliton 
{ki = —0.25, K,2 = 0) with small perturbation along the J = 1 
unstable eigenvector. Cross-sections of field intensities \En\'^ 
are plotted for the 0th (left column) and 1st (right column) 
harmonics at different propagation distances z. The initial 
soliton has a "donut" shape, which is consistent with the vor- 
tex charge /i = 1 in the first harmonic. As the instability 
evolves, it deforms the excitation towards usual bright spatial 
soliton bearing no vortex charge. Note, that the overall or- 
bital angular momentum is conserved and carried by rapidly 
diffracting radiative waves. 
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FIG. 6: (color online). The same as Fig. [5] but with pertur- 
bation along J — 2 unstable eigenvector. As the instability 
evolves, the soliton is transformed into a pair of spatial soli- 
tons, which are then pulled apart and moving in opposite 
directions. 





FIG. 7: (color online). Region of existence of coupled vortex 
solitons for the case of three fields: (a) asymmetric configu- 
ration (1-2 = -2,/_i = -l,/o = 0), [3-1 = 0.005, (3-2 = 0.02; 
(b) symmetric configuration (/_i = — l,/o = 0,/i = 1), 
(3i — [3-1 — 0.005. Open squares, filled circles and filled trian- 
gles correspond to numerically found instability thresholds for 
J = 1, J — 2 and J = 3 unstable perturbations, respectively. 
Shaded areas indicate regions of unstable solutions. 



by the condition = 1/2 is now K2 < i^2^\ where 



hereC^ = l/{/i2 + |^^_^^*_^ + ^2-m^3-mI'}- Fixing 
^1 — ^0 = 1, we consider two cases (M = 3 and M = 2): 
asymmetric {I-2 = — 2,/_i = — l,/o = 0) and symmetric 
(/_i = —1, Iq = 0, /i = 1). The former corresponds to the 
often encountered case with negligible ant i- Stokes side- 
bands, and the latter implies that the first Stokes and 
first anti-Stokes lines are excited. 

The existence boundary for the asymmetric case given 



= -]{—- 3ki - 2/3-1 - 13.2 
4 [ki 

In the symmetric case, the q = 1/2 condition implies 
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< K2 < K. 



(+) 



where 



P+ (1 - 4^2) 



4^2 



(25) 

with P± = ± Together with the conditions in 
Eqs. ([8|, the above constraints define the regions of the 
soliton existence, see Fig. ^ 

Stability analysis demonstrates that, similar to the 
two-component case with Iq = and /i = 1, the three- 
component solitons with l_2 = —2, l-i = —1, Iq = are 
stable inside a sufficiently wide domain in the {hZi,hZ2) 
plane and, in particular, in the proximity of the existence 
boundary given hy q = 1/2, i.e. in the high saturation 
regime. Close to the lower boundary of the existence do- 
main given by hZ2 = (^i + /5-2)/2 there are three types 
of instabilities with J = 1,2,3, see Fig. [7|^a). We note 
that the solution with the side-bands generated on the 
anti-Stokes side, i.e. the solution with /q = 0, /i = 1, 
/2 = 2, has the same stability properties as the solution 
discussed above. The symmetric case with n = —1, 0, +1 
is found to be unstable with respect to the J = 1 and 
J = 2 instabilities, with the former one persisting in the 
entire existence domain, see Fig. [7|^b). 



VII. MULTI-COMPONENT VORTEX 
SOLITONS AND SPATIO-TEMPORAL HELICAL 
BEAMS 

The above results show that if the vortex soliton con- 
tains a vortex free component, for example, at n = 0, 
and vortex carrying side-bands either only on the Stokes 
or only on the anti-Stokes sides, it can be stable within a 
broad range of parameters hii^2 ensuring that the satura- 
tion effects are sufficiently strong. Since in the frequency 
comb generation experiments with the off-resonant Ra- 
man gases the total number of excited harmonics can 
go to a few dozen [9l |1^, an important question to be 
addressed is whether the above stated principles of the 
vortex soliton stabilization can be extended onto multi- 
component cases. To address this problem we use numer- 
ical integration of Eqs. ([T]) with 11 coupled side-bands, 
initialized with the three-component vortex solitons de- 
scribed in the previous section. We consider two cases: 
(i) asymmetric case where excitation of the anti-Stokes 
lines is suppressed, and (ii) symmetric case with excita- 
tion of Stokes and anti-Stokes lines being equally proba- 
ble. In both cases we number the harmonics in a way that 
n = corresponds to the vortex- free component. Thus 
we take M = 11, A/" = and M = 6, TV = 5 in Eqs. Q 
for the asymmetric and symmetric cases, respectively. 

We monitor the evolution of the fields by plotting the 
total field intensity Itot = \Etot\'^ with Efot defined in 
Eq. ([2|. It has been demonstrated in [TT] that simul- 
taneous frequency and vortex combs lead to the helical 
structure of the total field intensity /tot, both in 




FIG. 8: (color online). Dynamics of 11 coupled fields initially 
excited with 3-component solitons. (a) Itot{x,y, z,t = 0) 
isointensity (x, 2;)-plot at 80% at maximum for asymmet- 
ric configuration. Fields n — —2,-1,0 are initialized with 
the soliton, ki = —0.25, Hi2 = 0.7 [stable for 3-component 
configuration, cf. Fig. |7|^a)]; (b) The same as (a) but for 
symmetric configuration, fields n = —1,0, 1 are excited with 
the soliton, k,i = —0.25, ^2 = [unstable for 3-component 
configuration, cf. Fig. [TJ^b)]. Isointensity plot is at 60% at 
maximum; (c) and (d) intensity distribution over harmonics 
after propagation distance z = 40 for the cases in (a) and (b) , 
respectively. 
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FIG. 9: (color online). Intensity distribution in the transverse 
plane of the 0th harmonic after propagation distance z — 
200 for the asymmetric (a) and symmetric (b) configurations 
in Fig. [S] Impact of the J = 3 instability in the case of 
symmetric configuration is clearly observed. 



and (x, z) subspaces. For the case of initial conditions 
where all the fields apart from the three pumps /c — 1, 
k and k -\- 1 are initially zero, the Itot can be crudely 
approximated with [see Appendix for details]: 

hotix, y, z, t) « |/f ) + fi°!,e-'* + /f+\e'*|' , (26) 

where (j) = t -\- AW — Kz^ Al = — h is the vortex 
charge step between the neighboring side-bands and K = 
uOrnodL/c — ^2- For any fixed t and z the total intensity 
distribution in the transverse plane is modulated in 
with the period defined by A/, and it rotates in both t 
and z^ forming a spatio-temporal helix. 

A/ = 1 corresponds to the single-strand helical struc- 



7 




200 




200 



FIG. 10: (color online). Stable double- (a) and triple-strand 
(b) helical beams formed by solitons with A/ = 2 and A/ = 3, 
respectively. All the parameters are the same as in Fig.pTa). 



ture of /tot, see Fig. [s] (a) and (b). Fig. [s] (a) shows 
the long distance evolution of the helix in the case of 
the asymmetric excitation, with all the side-bands gener- 
ated on the Stokes side, see the corresponding spectrum 
in Fig. [s] (c). The resulting helix in this case keeps its 
structure fixed over considerable propagation lengths. A 
similar numerical experiment for the symmetric excita- 
tion results in the helical soliton, which breaks up into 
filaments after the same propagation distance, cf. Figs. 
8(a) and (b). Note, however, that the total length in the 
simulations shown in Fig. [8] corresponds to a physical 
distance of order 20cm, which implies that one can speak 
about a quasi-stable propagation of the helix even in the 
case of the symmetric excitation of the Stokes and anti- 
Stokes side-bands. The z-period of the helix, 27r/K, is 
not a parameter of our numerical model, and it is only 
important when we are calculating Etot- Physically real- 
istic values of the adimensional period are of the order of 
1 (for a typical modulation frequency uJmod of the order of 
100 GHz [9 ), which makes the helical structure contain 
several hundred periods over the distance of 180 adimen- 
sional units required to see the instability. Therefore, to 
make the structure of the helices and the break-up pro- 
cess more obvious to the reader, we have fixed K ^ 0.1, 
when we have been producing the images of the helices 
in Figs. [8] and [To 



Providing the asymmetric excitation conditions and 
changing A/ to 2 and 3, we have also observed the for- 
mation of the stable double- and triple-strand helices. 



see Fig. 10 Note that the formation of similar multiple- 
strand helices has been reported in [25 , as a result of the 
linear superposition of the higher order Laguerre-Gauss 
modes. The helical soliton beams reported here are qual- 
itatively different from the so-called spiraling solitons or 
rotating soliton clusters |26l [271 EH EH EO] , which sustain 
their rotation due to the interaction between the individ- 
ual beams accompanied by the conservation of the an- 
gular momentum. In our case the helical evolution does 
not require the presence of more than one intensity lobe, 
as shown in Fig. [Sj and originates from the interaction 
of multiple frequency harmonics carrying progressively 
growing vortex charges. Most close known to us ana- 
logue of the spatio-temporal helices studied above have 
been reported in the context of the sine-Gordon equation 
and can be observed in a chain of coupled pendulums [31] . 



VIII. SUMMARY 

In this work we have reported existence conditions and 
have carried out linear stability analysis of the two and 
three component vortex solitons in an off-resonant Ra- 
man medium. We have found that, in the case where 
the vortex carrying Raman side bands are located either 
only on the Stokes or only on the anti-Stokes side of the 
vortex free component, the vortex solitons have a signifi- 
cant stability domain, corresponding to parameter values 
ensuring sufficient levels of nonlinearity saturation. We 
have also demonstrated that the same stabilization mech- 
anisms work in the case of many side-band, leading to the 
excitation of stable helical beams with single-, double-, 
and triple-strand topologies. 



Appendix 

An approximate expression for the z— evolution of the 
simultaneous frequency and vortex combs, excited with 
finite number of the side-bands, can be found if one ne- 
glects diffraction and dispersion. We replace Eq. Q with 
En{x^ z) ~ fni^)^^^"^^ and use the fact that under these 
approximations 



^(/n+l - Jn-l) 



(27) 



A solution to an initial value problem for Eqs. (27) can 
be expressed using the Bessel functions Jn{z). For an 
initial excitation with A^o adjacent side-bands: fn"^ 7^ 



for n = k^k 
given by 

fniz) 



1, 



A^o — 1, the resulting solution is 



k+No- 

E 

j=k 



^("-^■)/V„_,(29o^) , (28) 
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where qo = q{z = 0). The simplest case A^o = 2 has 
been considered in [3 [11] . Using the orthogonaUty of the 
Bessel functions: Jn+pJn+g = ^p,qi it is easy to show 
that q{z) = qo and thus Eq. (28) satisfies Eqs. ([27| for 
ah z. Substituting the solution (28) into Eq. ([2|, we find 



the approximate expression for the total field: 



ino\ X 



Etot - exp(i^o) {^^P [' 

n 



where = ^o^ + ^ot/^mod — Kqz, <j) = /\16 + t — 
Kq — (JqL/c — /^i, K — iJrnodL/c — K,2' Usiug a kuowu 
identity, Jn{x) ex.p{ina) = exp[zx sin(a)], we derive 



/c+A/'o-l 

Etot ^ exp[i(l)o^i2qozcos{(l))] ^ exp[ij^], (30) 

j=k 



(29) which is the expression used in Eq. (26). 
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